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Abstract 

The superheating that usually occurs when a solid is melted by vol¬ 
umetric heating can produce irregular solid-liquid interfaces. Such inter¬ 
faces can be visnalised in ice, where they are sometimes known as Tyndall 
stars. This paper describes some of the experimental observations of Tyn¬ 
dall stars and a mathematical model for the early stages of their evolntion. 
The modelling is complicated by the strong crystalline anisotropy, which 
results in an anisotropic kinetic nndercooling at the interface; it leads to 
an interesting class of free boundary problems that treat the melt region 
as inhnitesimally thin. 


1 Introduction 

When a single crystal of pure, transparent ice is irradiated, the partial absorp¬ 
tion of transmitted radiation volumetrically heats the crystal, leading to internal 
melting and the formation of small volumes of liquid. Remarkably, these vol¬ 
umes of water often take on shapes that resemble six-fold symmetric flowers, 
stars, or snowflakes, as first documented by Tyndall m- The internal melt 
figures that Tyndall observed now bear his name and are often referred to as 
Tyndall stars, Tyndall figures, or liquid snowflakes. An examples of such can 
be found in Fig. 

Tyndall stars are predominantly found in very pure crystals of irradiated ice. 
The lack of impurities and microscopic defects in such crystals limits the onset 
of liquid nuclei and prevents the ice from simply melting away as it continually 
absorbs radiation. Instead, the ice becomes superheated, whereby its tempera¬ 
ture exceeds the equilibrium melting temperature. It is this superheating that, 
through an interfacial instability, is suspected of giving rise to the complex mor¬ 
phologies that are characteristic of Tyndall stars. The six-fold symmetry that 
is apparent in Fig. is inherited from the anisotropy of the ice crystal, which 
will be discussed in detail below. 


1 



Figure 1: An example of a Tyndall star that has been created by irradiating 
a pure crystal of ice with light from an overhead projector. The bright circle 
within the star is a vapour bubble that emerges due to the density difference 
between water and ice. The viewing plane corresponds to the basal plane of 
the melting ice crystal with the c axis pointing orthogonally into and out of 
the page. This image was created by the authors at the FoaLab in Oxford; for 
additional details, see Harvey M- 


From a scientific viewpoint, Tyndall stars offer a convenient route for study¬ 
ing the dynamics of phase change and moving interfaces because both the solid 
and liquid phases are transparent. Thus, in principle, these phases can be 
observed in real time with visible light. Understanding of Tyndall stars may 
also have industrial implications in, for example, resistance welding, whereby a 
metal is volumetrically heated by passing an electrical current through it min]. 
This leads to a superheated solid and the formation of small inclusions of liquid 
metal. Due to the opacity of the metal, these inclusions cannot be seen in real 
time and are often detected after the welding operation is over. 

The evolution of Tyndall stars has been studied experimentally by Nakaya 

, who found that the melts begin as cylindrical discs of water with thicknesses 
that are much smaller than their radii. This thin aspect ratio is maintained 
during the evolution of a Tyndall star, with growth in the radial direction being 
much faster than in the axial direction. As the cylindrical disc increases in 
size, the circular interface can become unstable, leading to the emergence of a 
high-wavenumber fingering pattern. In cases where the radiation intensity was 
sufficiently high, further growth of the instability resulted in the formation six 
large symmetric dendrites. In addition, Nakaya reported that the Tyndall stars 
in a given ice crystal always have the same orientation. Further experiments by 
Takeya [29] were able to provide quantitative data for the radial and axial growth 
of Tyndall stars. Over the duration of a couple of minutes, the radius increased 
to roughly 1.5 mm while the thickness grew linearly with time to about 0.3 
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mm. In some cases, however, the axial growth of the melt was only temporary 
and eventually it stopped altogether. Interestingly, Takeya reported that an 
interfacial instability only occurs when the axial growth persists; in cases where 
the axial growth terminates, the melt remains cylindrical^ This observation 
is perhaps linked to those made by Mae [52], who found that Tyndall stars 
retain their initial cylindrical shape unless they grow beyond a critical thickness 
of 10 /rm. Experimental |2H] and theoretical m studies of solidification in 
supercooled liquids, a situation that closely parallels melting into a superheated 
solid, have also shown that a critical thickness must be surpassed in order for a 
morphological instability to occur at the solid-liquid interface. 

The anisotropic growth of a Tyndall star is closely related to the geomet¬ 
ric configuration of the melting ice crystal. Roughly speaking, the crystalline 
structure of ice can be imagined as a collection of adjacent hexagonal prisms; 
see Fig. The hexagonal faces of the prisms form the so-called basal planes 
of the crystal and the direction that is normal to these planes defines the c 
axis. The radial growth of Tyndall stars occurs within the basal planes while 
the axial growth is aligned with the c axis, therefore giving different Tyndall 
stars the same orientation within an ice crystal. The molecularly smooth basal 
planes melt at a much slower rate than the molecularly rough prism planes. As 
discussed in the context of solidification |5] , the accretion of material normal to 
a molecularly smooth surface, such as a basal plane, occurs via an energetically 
activated process, whereas there is no nucleation barrier at a molecularly rough 
surface. The fast-melting prism planes dominate the shape of the Tyndall figure 
|25j and are responsible for the disparity between its axial and radial dimensions. 

The mathematical study of problems involving phase change is now a clas¬ 
sical subject for which there is extensive literature. Davis nni gives a compre¬ 
hensive treatment of the mathematical theory of solidihcation, starting from 
the classical Stefan problem. Hu & Argyropoulos [13] provide an overview of 
modelling and computational techniques that are relevant to solidification and 
melting problems. The fluid mechanics of solidification are reviewed in detail by 
Huppert m- Coriell et al. 015] examine the occurrence of multiple similarity 
solutions, as well as their selection mechanisms, in models of solidification and 
melting. The application of phase-field models to solidification problems has 
been discussed by Boettinger et al. [3]. 

Mathematical models of phase change that account for the anisotropic na¬ 
ture of the solid have been largely confined to the case of solidihcation and 
crystallisation. Wettlaufer et al. [6ll23l|3ll[33] examined two-dimensional crys¬ 
tallisation within the basal plane by considering an interfacial velocity that 
depends on the angle between the free boundary and a certain hxed direction. 
A suitable angular dependence was found to give rise to the six-fold symmetry 

^The axial growth ceased for cases of low superheating with there being sufficient heat to 
melt only a small part of the ice. There could be only limited scope for instability in such 
situations. 
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Figure 2: A schematic diagram of an ice crystal, which is composed of arrays 
of hexagonal prisms. Shaded hexagonal faces form the molecularly smooth 
basal planes of the ice crystal and the unshaded rectangular faces correspond to 
molecularly rough prism planes. The c axis of the ice crystal is orthogonal to 
the basal planes. The shaded circles give the approximate positions of oxygen 
atoms. The rate of melting is much higher at prism planes than basal planes, 
resulting in Tyndall stars that are relatively thin in directions along the c axis. 


that is characteristic of snowflakes. It is important to emphasise here that in 
the studies of Wettlaufer et ai, it is assumed that growth of the crystal is in the 
geometric limit, whereby the interface velocity is only a function of the shape 
and position of the interface. In particular, the velocity of the interface does 
not depend on field variables that are affected by its motion. This is in contrast 
to non-geometric growth models, which account for long-range diffusion of field 
variables and their coupling to the interfacial velocity. In geometric models, the 
crystalline anisotropy enters directly through the interface velocity. However, 
in non-geometric models, anisotropy enters through physical parameters related 
to the interface, such as surface energy or the coefficient of kinetic undercool¬ 
ing, the latter of which connects the temperature and velocity at the interface. 
Anisotropic solidification outside of the geometric limit has been investigated by 
a number of authors. Uehara & Sekerka [32] studied the formation of facets due 
to strong anisotropy in the kinetic coefficient using a phase-fiefd modef. Par¬ 
ticular attention was paid to determining the relationship between the shape of 
the emerging crystal and the mathematical properties of the anisotropic kinetic 
coefficient. Yokoyama & Kuroda [35] employed the boundary-element method 
to study the hexagonal morphologies of snow crystals predicted by a model 
with an anisotropic kinetic coefficient. Yokoyama & Sekerka |3bj explored the 
combined effects of anisotropic kinetic undercooling and surface energy. Using 
numerical and asymptotic methods, they investigated the suppression of corner 
formation between adjacent facets. 

Considerable attention has focused for many years on the stability of the free 
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boundary in phase-change models. Linear stability analyses of models which 
treat the phase interface as infinitesimally thin, such as in the pioneering study 
by Mullins & Sekerka [26] or in Hele-Shaw and Muskat problems, indicate that 
a morphological instability can arise when a melting boundary is driven by 
heat flow from a superheated solid region |20j . In fact, without a regularising 
mechanism such as surface energy or kinetic undercooling, the system is severely 
unstable and the model becomes ill posed in the sense that disturbances with ar¬ 
bitrarily large wavenumbers will grow arbitrarily fast in time. Such ill-posedness 
can also be avoided by replacing the sharp, infinitesimally thin interface with a 
diffuse mushy region consisting of two co-existing phases niinj. The theory of 
mushy regions in volumetrically heated solids has been developed by Lacey et 
al. naEani], who treated the mush as a collection of small liquid inclusions 
that grow within the solid. In these papers, the growth of the inclusions is 
modelled using classical Stefan problems that account for surface-energy effects 
and interfacial curvature, kinetic undercooling, and/or composition in the case 
of alloys. The main purpose of those studies was to use homogenisation to build 
an averaged model for the mushy region. 

A sharp-interface model of Tyndall stars has been formulated and studied 
by Hennessy m- The focus here was on two-dimensional evolution within the 
basal plane. The morphology of the solid-liquid interface was studied using a 
combination of linear stability theory and numerical simulations. Growth along 
the c axis was not considered and thus it was not possible to explore how this 
may influence the stability of the ice-water interface. 

In this paper, we consider the three-dimensional evolution of a Tyndall star 
or, perhaps more accurately, a Tyndall hgure, as we mostly discuss the ear¬ 
lier growth rather than the later, star-like stage. Particular attention is paid 
to capturing the anisotropic growth along the radial and axial directions. Our 
description of the problem is based on the classical Stefan model but the inclu¬ 
sion of volumetric heating and anisotropic kinetic undercooling makes it non¬ 
standard. An asymptotic analysis that exploits the axial and radial length-scale 
separation is used to reduce the three-dimensional problem to a co-dimension-2 
free boundary problem whereby the melt is collapsed into a planar surface with 
infinitesimal thickness. A local stability analysis of the reduced model is carried 
out as a first step towards the study of the onset of fingering patterns at the ice- 
water interface. An attempt is made to compare our theoretical results to the 
experimental observations of Takeya |29j : however, this is not straightforward 
due to a lack of knowledge of key quantities controlling the anisotropic growth. 
We then propose future experiments that could produce novel quantitative in¬ 
sights into the growth kinetics. 

In the next section, we present a mathematical model for a growing Tyndall 
figure based on laboratory experiments. In Sec. we carry out an asymptotic 
analysis of this model that captures the anisotropic growth of the melt and 
investigates the stability of the ice-water interface. We discuss our results and 
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Figure 3: We study the growth of a Tyndall figure (depicted by the shaded 
region) in superheated irradiated ice. We use and to denote regions of 
space occupied by liquid water and solid ice, respectively. Here, t represents 
time. The ice-water interface is denoted by T(t) and has a normal vector n and 
normal component of velocity v. The z axis is parallel to the c axis of the ice 
crystal and r = i® ^ radial coordinate that lies within the basal 

plane. The angle between the c axis and the normal vector n is given by ip. 


conclude the paper in Sec. 


2 Mathematical Model 

2.1 The Physical Problem 

We suppose that a single crystal of ice held at its melting temperature is il¬ 
luminated at time t = 0. The direction of the incident light is taken to be 
parallel to the c-axis of the crystal; see Fig. We assume that a rapid nu- 
cleation process occurs within the ice upon exposure to light, leading to the 
creation of a single spherical melt figure. Continued absorption of radiation by 
both the ice and the water will drive the melting at the interface, which we aim 
to describe mathematically. Our model of this physical scenario is based on 
equations governing the temperatures in the liquid and solid phases, taking into 
account thermal diffusion and volumetric heat generation due to absorption of 
radiation. The solid-liquid interface is assumed to be sharp and, therefore, we 
impose appropriate boundary conditions on it. 

The field equation for the temperature Tj of phase j is given by 

dr 

+qj, X e %(t) (1) 

where t is time, position is a; = (x, y, z), and is the region of space occupied 
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by phase j. We let j = I and j = s for the liquid water and solid ice phases, 
respectively. We assume that the z axis and the {x, y) plane are aligned with the 
c axis and basal planes of the ice crystal, respectively. The values of the material 
constants, namely the densities, pj, specific heat capacities, Cpj, and thermal 
conductivities, kj, differ between the two phases. Although the difference in 
density between the phases is significant enough to give rise to a vapour bubble 
inside the Tyndall figure, as shown in Fig. their relative difference is small 
and we take the densities of the two phases to be the same and equal to p, that 
is. Pi = Ps = p. The rates of volumetric heating, are given by the product 
of an absorption coefficient, p,j, and the local intensity of incident light upon 
the medium, I. With a sufficiently small piece of ice (or absorption coefficient), 
I can be regarded as constant, making qj constant in each phase. We shall 
generally assume that the initial temperatures coincide with the equilibrium 
melting temperature Tq at t = 0, with a spherical Tyndall figure of radius a 
nucleating at the same instant. However, if significant body heating occurs 
before nucleation, the initial temperatures will be much greater than Tq. This 
situation is discussed in Appendix [B} 


At the evolving interface T = T{t) between ice and water, we have the usual 
Stefan condition 


Lpv = 


dT^ 

dn 


e r(t). 


( 2 ) 


where L is the latent heat of fusion, assumed constant; v is the normal velocity, 
measured towards the ice; djdn is the normal derivative, again in the direction 
into the ice; and [-J® denotes the change in a quantity across the interface, going 
from liquid water to solid ice, see Fig. 


We also assume that the normal velocity of the interface is proportional to 
the local amount of superheating m- To account for the different melting rates 
of the basal and prism planes, we take the constant of proportionality to be a 
function of the orientation of the interface. Thus, we impose a kinetic condition, 
equivalent to anisotropic kinetic undercooling in solidihcation |32[ 1351136j , given 

by 


v = Kf{ij){Ti-To), xeT{t), (3) 

where T/ is the temperature at the interface, 

Tj=Ts = Tu x&T{t)- (4) 

K is a constant; f is a dimensionless function, which we refer to as the anisotropy 
function; and ip is the angle between normal vector at the free surface and the 
c axis, as measured relative to the positive x axis, see Fig. Contributions to 
<§ from the surface energy are not included, which we justify by assuming that 
after the rapid nucleation phase, the radius of the melt is much larger than the 
capillary length given by leap = il/pL){To/AT), where 7 is the surface energy 
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of an ice-water interface and AT is the local amount of superheating. Takeya 
[221 measured superheatings on the order of 0.1 K in his experiments that use 
photographic bulbs as the light source, giving a capillary length of 270 nm; thus, 
neglecting surface energy seems reasonable given that Tyndall figures typically 
have length scales on the order of hundreds of microns up to millimetres. An 
important consequence of neglecting surface energy is that our model will not 
capture the evolution of the melt into a Wulff shape, which is the equilibrium 
shape arising from the minimisation of surface energy under constant-volume 
conditions [10]. However, based on the phase-field simulations by Uehara & Sek- 
era , we might expect the melt to grow into its “kinetic Wulff shape”, which, 
in essence, describes the asymptotic shape that the interface would approach if 
it were to evolve solely due to anisotropic undercooling under isothermal con¬ 
ditions, so that the normal velocity depends only upon the orientation of the 
interface [32l[35l|36] (also see, below, Sub-sec 3.1 and Sub-sec 3.2). 


2.2 The Anisotropy Function 

The anisotropy function / is used to model the orientation dependence of the 
interfacial velocity arising from the crystalline structure of the ice. We assume 
that the value of / is close to one when the velocity is parallel to the prism planes 
of the ice crystal and small when the velocity is parallel to the basal planes. 
Mathematically, this corresponds to / ~ 1 when ip = ±7r/2, and / ~ e ^ 1 
when Ip = Q, ±7r, respectively. In physical terms, the parameter e can be thought 
of as the ratio of the melting velocity of basal planes to prism planes for a fixed 
superheating Tj — Tq > 0. Experimentally determining a functional form for / 
is possible by measuring the kinetic Wulff shape. However, acquiring the kinetic 
Wulff shape is difficult in practice and, consequently, there is often uncertainty in 
the form of /. Therefore, our analysis will rely on phenomenological expressions 
for the anisotropy function. More specifically, we will consider in detail the 
function 

f{iP) = ie^ + sin^iP)y^, (5) 


which is expected to produce smooth interfaces based on its corresponding ki¬ 
netic Wulff shape. In two-dimensions, the kinetic Wulff shape is determined by 
the convex region containing the origin traced out by the parametric curves 


X = f {ip) cosip + f (ip) sin Ip, (6a) 

z = f {ip) sin Ip — f {ip) cosip. (6b) 

The anisotropy function ([^ is shown along with its corresponding Wulff shape 
in Fig. Additionally, we will present the key results that are obtained when 

f {ip) = e + sin'^ Ip, (7) 


and 


/(V') 


e 

1 -I- e — sin^ Ip 
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Figure 4: Top (a): we consider three different anisotropy functions / that char¬ 
acterise the dependence of the melting rate on the orientation of the solid-liquid 
interface: /(i/^) = (e^ -|-sin^^)^/^ (solid), /(^/>) = e -I- sin^ %j) (dash-dotted), and 
= e/{l + e — sin^ ijj) (dashed). Here, ip measures the angle between the c 
axis and the vector normal to the interface; see Fig. Bottom (b): the corre¬ 
sponding kinetic Wulff shapes associated with the three anisotropy functions /, 
which represent the long-term shape the melt would acquire under isothermal 
conditions and growth due purely to anisotropic kinetic undercooling. See text 
for further details. In both panels we have set e = 0.1. 


Since sin^ ip can be written in terms of sin 2ip, the anisotropy function Q is 
similar to many of those found in the literature [32]. The anisotropy function 
in Q has sharp maxima ad, ip = ±7r/2 (see Fig. (a)), making it comparable 
to theoretical expressions for / that have been derived from models of surface 
diffusion 13135]. Figure |4] shows that the anisotropy functions 0 and ^ lead 
to the formation of corners in the kinetic Wulff shape. Surface energy is likely 
to become important on these small scales and may lead to a smoothing of the 
corners. Capturing such dynamics is beyond the scope of our current model, 
however. 

2.3 Parameter Values 

The configuration that we study here is based on experiments involving ice-water 
systems carried out in Oxford. Light from an overhead projector was used to 
irradiate a pure ice crystal. A list of parameter values corresponding to these 
experiments is given in Table[^ Although light from the overhead projector will 
have a broad spectrum, ice and water are particularly strong absorbers of infra- 
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Table 1: Parameter values for ice-water systems heated by light from an over¬ 
head projector. These are based on experiments carried out in Oxford. The 
absorption coefficients are for monochromatic infra-red radiation with a wave¬ 
length of 980 nm. The intensity of radiation is estimated from the power of the 
bulb and distance to the sample, further details are given in the text. 

p 1000 kg/m^ 

Cps 2050 J/(kg K) 

Cpi 4181 J/(kg K) 

ks 2 W/(m K) 

ki 0.6 W/(m K) 

L 3.33 X lO'^ J/kg 

/is 15.3 1/m 

Pi 43.6 1/m 

lo 300 W/m2 

7 0.033 J/m2 

To 273 K 


red radiation. Therefore, the absorption coefficients in Table are based on 
monochromatic infra-red radiation with a wavelength of 980 nm. The intensity 
of radiation has been calculated from the bulb power and distance to the sample 
by assuming spherical emission; the complete details can be found in Hennessy 

m- 

Determining values for the parameters K and e is a challenging experimen¬ 
tal task. Using arguments from statistical mechanics, it is possible to write the 
velocities of the planes [34], as well the coefficient K in (§ m , in terms of 
elementary quantities such as molecular distance and activation energy. How¬ 
ever, these expressions introduce additional unknown parameters into the model, 
making them of little practical use. The combined uncertainty in the values for 
K and e, as well as in the functional form of the anisotropy function /, will 
make carrying out a quantitative comparison of our results with experimental 
data difficult. That being said, qualitative comparisons are still possible, and 
the analysis can be used as a tool for ruling out anisotropy functions. 

2.4 Non-dimensionalisation 

The model is non-dimensionalised by introducing suitable scales for time, dis¬ 
tance, and temperature. The time variable t is written in terms of the time 
scale of thermal diffusion in ice, jus-, where Ks = ks/{pcps) is the thermal 
diffusivity of ice and £ is a characteristic length scale defined below. The tem¬ 
perature scale is set by the amount of superheating in the ice caused by vol¬ 
umetric heating, giving AT = /kg- Finally, the length scale ^ is chosen 
to balance terms in the kinetic condition ([^, implying that significant growth 
parallel to the basal planes occurs on 0(1) (dimensionless) time scales. This 
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gives = Kk1/{qspcps). Using these scales, we write t = {£'^Iks)t, x = iX, 
and Tj = Tq + {AT)9j. The non-dimensional field equations can be written as 

X G r!,(T), (9a) 

X G ni{T), (9b) 

where Cp = Cpifcps and k = ki/ks are ratios of specific heat capacities and 
thermal conductivities, respectively. The ratio of volumetric heating, q = qi/qs, 
can be written in terms of the absorption coefficients via q = pi/ps- Initial 
conditions for the temperatures are given by 0s = 0; = 0 at r = 0. 

At the ice-water interface, the Stefan and kinetic conditions, along with the 
continuity of temperature, are given by 


, (do, .de,\ 

, XGr(r), 

(9c) 

V = 01 fW, 

X G r(r), 

(9d) 

II 

Cfi 

II 

X G r(T), 

(9e) 


+ 1 , 

OT 

Cp^ = ky^9i + q, 


respectively, where j5 = L/{cpsAT) is the Stefan number. The initial ice-water 
interface is taken to be the sphere with dimensionless radius a = aj£ given by 
|X|=a. 

Far from a growing liquid inclusion, d9s/dT ~ 1, so that we have 

0s ~T, |X|^C30. (9f) 

Note that (|^ requires that the Tyndall figure and associated length scales be 
small compared with the region subject to the body heating. 

Using the parameter values in Table we find that k ~ 0.3, Cp ~ 2, and 
g ~ 3, all of which can be treated as 0(1) in size. Due to uncertainty in 
the value of the parameter K, it is difficult to estimate the length scale £, the 
characteristic temperature rise AT, and the Stefan number f3. Using instead the 
measured value of AT ^ 0.1 K from Takeya [29], the Stefan number is given by 
P ~ 10^. The length scale can be estimated from I = (ATfcs/gs)^^^ 6-7 mm 
and the time scale from €^/ks ~ 46 s, which seem slightly large but reasonable. 

The proceeding analysis will focus on the distinguished limit whereby e = 
0{P~^). This regime is considered so that we can examine the interplay of the 
kinetic anisotropic effects; whether or not this balance occurs in practice depends 
upon the size of the rate of the volumetric heating. Thus, we write = be 
where b = 0(1). Furthermore, it will be assumed that the (dimensionless) radius 
of the initial melt, a, satisfies a <?; e. In dimensional terms, this inequality 
means that the initial radius should be less than one micron, which is close to 
the limit where surface-energy effects become important. This upper bound on 
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the initial size of the radius, along with the anisotropic kinetic condition (9e), 
ensures that the spherical Tyndall figure will first grow into a thin disc of melt 
with radius that is much greater than its thickness, which is consistent with 
experimental observations |29j . 


3 Analysis 


The analysis begins in Section |3.1| with an examination of the small-time be¬ 
haviour for r = In dimensional terms, the small-time regime corre¬ 

sponds to times given by t ~ {ajjK s). Taking the dimensional radius of 
the initial melt to be of the order of one micron, we find that t ~ 0.5 seconds. In 
this regime, the volumetric heating and the kinetic condition drive the melt into 
a thin shape with dimensions along the c axis that are much smaller than those 
in the basal plane. In Section 3.2 we consider the dynamics when t = 0(1), 
corresponding to t ^ 50 seconds. By exploiting the separation of length scales 
that arises from the initial growth, a simplified model can be derived. Using this 
model, the linear stability of the ice-water interface is examined in Section |3.3[ 
Our analysis will first focus on the dynamics that occur when the anisotropy 
function ® is used. We will then consider additional anisotropy functions in 
Section 13.^ 


3.1 Early Time 

The analysis of the early-time behaviour proceeds by letting r = 6j = 

where a ^ e 1. We then consider the temperature field near and away 
from the melt, and connect the solutions in the two regions using asymptotic 
matching. 


In the region of solid away from the melt, i.e., for X ~ 0(1), the leading- 
order problem in a is straightforward to solve and it gives 9s{X,t) = f. To 
resolve the temperatures near the melt, we let X = aX. The leading-order 
problem in a in this inner region is given by 

\ 7 ^ 9 , = 0, X e U«(f), (10) 

v¥i = o, xeCiiif), (11) 


with the following conditions at the solid-liquid interface: 


dn dn’ 


xef(f). 


( 12 ) 

(13) 


By asymptotically matching the temperatures in the solid, we also have that 
—)• T as \x I —>■ oo. The solutions for the temperature fields are given by 
9i = ds = f. The motion of the interface, therefore, satisfies the equation 


V = rfiip). 


(14) 
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To make further progress, we suppose that the rescaled positions of the ice- 
water interface are given by the zero level set of a function F, defined by 


F = s{X{f))-f^/2 = 0, 


(15) 


where, s a function that is to be determined. The initial shape of the interface 
is encoded in the function s; we require that s(.X'(0)) = 0 on the sphere |.X’| = 1 
when T = 0. We emphasise here that s also plays the role of a time variable; 
from (15) we see that s = f^/2. In this formulation, the normal velocity at the 
interface can be written as u = f/lVsl and, therefore, the kinetic condition (14) 
becomes 


\S/s\fW = 1 . 


(16) 


Closing the problem requires writing the angle ip in terms of the function s. 
For clarity, we now consider the two-dimensional problem by writing X = 
{X,0,Z). In this case, simple trigonometry shows that the angle ip satisfies 
sinip = sjf/(s^ -b where Sjy = ds/dX and = dsjdZ. By writing 


f{ip) = f {sinip), the kinetic equation (16) becomes 


(s^ + s| 


)‘'V( 




= 1 . 


(17) 


To see how the melt region evolves, we now focus on the anisotropy function 
given by ([^. In this case, the problem for s is 

(4 + 4 )'/= [,= + 4(4 + 4 )-.] = 1 , ( 18 ) 

subject to the condition s{Xo,Zq) = 0 on the circle X§ -|- Zq = 1 at time 
f = 0. The solution to this problem can be found using Charpit’s equations, 
as detailed in Appendix]^ In essence, Charpit’s equations are a generalisation 
of the method of characteristics for nonlinear first-order hyperbolic problems. 
We proceed by parametrising the initial data according to Xo((p) = (Aq, Zq) = 
(cost/?, sint/?), s(.X'o(t/3)) = 0, with ip € [0, 27r). Upon applying the method, 
solution can be written implicitly and parametrically as 


A = 


1 + 


s(I -I- e^) 


(e^ -I- cos^ 


cos ip, 


Z = 


1-b 


se^ 


(e^ -I- cos^ 


sint/?, (19) 


Thus, for a given value of s, which can be written in terms of time via s = 12, 

these curves trace out the instantaneous positions of the solid-liquid interface as 
ip is varied from 0 to 27r. Figure]^ shows the interface profiles predicted by ( |I^ 
at various times when e = 0.1. The initially spherical melt first grows primarily 
in the radial direction, keeping its thickness in the axial direction constant (top 
panel). By the time the axial growth becomes appreciable, the radius of the 
melt has grown a substantial amount, resulting in a liquid region with a small 
aspect ratio. 
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Figure 5: The early-time evolution of a spherical melt figure when the anisotropy 
function is given by ([^ when e = 0.1. These curves are given by the solution 
in (19). The arrows indicates the direction of time. The top panels shows the 
solid-liquid interface at equally spaced values of s given by s = 0, 0.33, 0.66, and 
1, corresponding to rescaled dimensionless times given by f = = 0, 0.82, 

1.15, and 1.41, respectively. Similarly, the bottom panel shows the interface for 
values of s given by s = 0, 1.42, 2.86, 4.29, 5.71, 7.14, and 8.57, corresponding 
to r = 0, 1.69, 2.39, 2.93, 3.38, 3.78, and 4.14. The interface remains smooth 
for all time and evolves into the kinetic Wulff shape shown in Fig. 
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To aid in the physical interpretation of (19), we revert to the original non- 
dimensionalisation by writing s = s/a, X = X ja, f = to obtain 


X = 


s(l + e^) 


(e^ + cos^ 


cosip, Z = 


se 


(e^ + cos^ 


sin(^, (20) 


where s = In the very early stages of development, so that s is of order 

a, then for parts of the interface given by | cos(/5| ^ e. 


Z ^ a sin p, X ~ a + 


cos p\ 


cos p ~ ±s + a cos p, 


while for | cost/?] = 0(e), say p = ±7r/2 =F with ip = 0(e), 

zL-ips 


Z ~ ±a, X ~ ± a ± 


■ 




(e2 +^2)1/2- 


( 21 ) 


( 22 ) 


Thus the interface takes the form, approximately, of two circular arcs, each of 
radius a and centred on (X,Z) = (±s,0), linked by horizontal lines. 


In the later stages, s » a/e, 


and 




X2 


(e^ + cos^ 


sin p, X 


s(l + e2) 


■ cos p, 


Z^ 
H— 


(e^ + cos^ 

[(1 + e^) cos^ p + e^ sin^ p] = s^. 


(23) 


(24) 


1 _|_ g2 g2 g2 _|_ (,Qg2 ^ 

so the interface is then approximately elliptical. The longer-term interface pro¬ 


file, defined by the large-time limit of the small-time model, and given by (23) 


for this choice of /, is, in fact, equivalent to the corresponding kinetic Wulff 
shape that can be computed from Note that the half thickness of the melt, 
given by the maximum value of Z, grows in time as Z(p = 7r/2) = eT^/2. 
The maximum value of X, corresponding to the rim of the melt, grows as 
X((/j = 0) ^ t^/2. We see that for s ^ a, the influence of the initial interface 
has been lost. 


3.2 Order-One Time 

We now consider the dynamics that occur on 0(1) time scales. The initial 
condition in this time regime takes the form of a matching requirement, as r —^ 


0, with the fully developed early-time shape given by (24). The analysis in two 


and three dimensions is sufficiently similar for us to proceed directly to problems 
with axial symmetry. Thus, we define a radial coordinate R — (X^ -|- 
We also assume the system remains symmetric about the Z = Q plane and, 
thus, we only consider the problem in the upper-half space given by Z > 0. The 
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Figure 6: A schematic diagram showing the three asymptotic regions in the 
T = 0(1) problem. By constructing local solutions in regions (i) and (ii), it 
is possible to derive effective boundary conditions that lead to a self-contained 
problem in region (iii) by asymptotic matching. 


position of the solid-liquid interface is written ss Z = h{R, r); the corresponding 
position of the rim is i? = S{t) so that = 0. The angle ^ appearing 

in the anisotropy function / satisfies 


. , dh/dR 

“ (l + (a/i/9A)2)i/2- 


(25) 


From matching into the early-time regime and using (|24[), we expect that 


/i(i?) ~ e(S'^ - i?2)i/2 


(26) 


as T ~ 0. 

In principle, the dynamics in the 0(1) time regime can be studied by solving 
§ directly. However, the thin aspect ratio of the melt, with Z ~ 0(e) and 
X,Y ~ 0(1), motivates seeking a solution via matched asymptotic expansions, 
and this is the approach we take. There are three distinct regions that need 
to be considered: (i) near the melt but away from the rim, (ii) near the melt 
and near the rim, and (iii) away from the melt. A schematic diagram of these 
regions is shown in Fig. Our approach is to obtain local solutions in regions 
(i) and (ii) which can then be used to derive effective boundary conditions for 
the problem in region (iii) by asymptotic matching. 

3.2.1 Analysis near the melt and away from the rim 

In region (i) near the melt but away from the rim, R <C S{t), we rescale the 
axial coordinate according to ^ = eZ. In addition, the position of the interface 
is written as h{R, r) = eh[R, r) and the temperatures in this region are denoted 
by 0j, j = Z,s. Under this scaling, the anisotropy function ([^ can be written 
as /('0) ~ e[l -I- {dh/dR^Y^^. The governing equations in this region are given 
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by 




dr RdR\ dR 


e^c, 


^0; 

dr 


e^d_ 

lilm 


R 


dQi 

OR 


J dZ^ 

Z > h{R, r), 

(27a) 

f d^Qi 2» 

k ~ + 6 Q, 

9Z2 

Z < h{R, t). 

(27b) 


The boundary conditions on the solid-liquid interface are 

,90, dh 


OT 


OT 


90; _ 2<^_ 

dZ ^ dR dR 


. / dh 

^ ( dR 


90; dh 
iWl 

2 -I 1/2 


90, 

dZ 


— , Z = h{R,T), (27c) 


'-Mi 


dR dR' 
2 -1 1/2 


Z = h{R,T), (27d) 


where 0/ = 0,(7?, h(i?, t), r) = 0;(i?, h(i?, r), r). The symmetry about Z = 0 
implies that dQi/dZ = 0 at Z = 0. The relevant matching conditions for the 
temperature in the solid as Z —>■ oo will be discussed below. 


The solution to this problem is now expanded as 

0, = ef + e©^ + Oie^) 
h = h^°'> + eh^^'> + 0{€^). 


(28a) 

(28b) 


Assuming that e^q = O(e^), the 0(1) solution for the temperature is straight¬ 
forward to obtain and is given by 


0™ {R,Z,T) = e<f\R,Z,T) = 0f (E, r). 


(29) 


The matching condition for this problem is given by ©1°^ (i?, Z, t) = 9s{R, 0, r) 
as .Z —>■ oo. From (29l, we can deduce that Q^P{R,t) — 0,(7?,0,t). Therefore, 


the 0(1) part of the kinetic equation (27d) becomes 
dh^°^ 


dr 


= 0,(7?,O,r) 


1 


9h0) 

97? 


1/2 


(30) 


Proceeding to the 0(e) problem, we find that the temperatures are deter¬ 
mined from bulk equations 


920 b 

- J— = 0, 

9Z2 


and must satisfy the Stefan condition 


,_,9h0) ~90) 

0 - = —k 




9t 


dZ 


+ 


dZ 


(31) 


(32) 
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By exploiting the symmetry of the problem about the Z axis, we find that the 
temperature in the liquid, must be constant in space. Asymptotically 

matching the derivatives of the solid temperature in regions (i) and (iii) gives 
the relation 


dZ dz 


(33) 


as 2’ —)■ oo and z 


0. Using ( |33[ ) in the Stefan condition (321 yields 
89, 


dr 


= b 


dz 


z = 0. 


(34) 


We emphasise here that (30) and (34) can be treated as boundary conditions 
for the problem in region (iii) away from the melt. 


3.2.2 Analysis near the melt and near the rim 


The next step is to consider the local dyiramics near the rim in order to derive an 
equation describing its motion. We switch to a travelling-wave coordinate given 
hy R= {R — S{T))/e^ and let Z = e^Z. These scales have been chosen in order 


to balance both sides of the initial interface profile given in (26). The position 
of the solid-liquid interface is written as h{R, r) = e^h{R) and the temperatures 
are denoted by Qj for j = I, s. Upon using this scaling in ([^, the leading-order 
problem in e is given by 


920 , 920 , 


9i?2 9Z2 


= 0, {R,Z) G j = s,l. 


(35) 


The Stefan condition reduces to the continuity of thermal flux across the inter¬ 
face: 


= Z = h{R). R<0. 

dZ dR dR \dZ OR dR ^ ’ 


The leading-order kinetic equation reads 

dh 

dr dR ~ ^ 


dR 


Z = h{R), R<0, 


(36) 


(37) 


where 0/ = Qs{R,h{T),T) = 0i(J?,/i(r), r). Since the thickness of the melt 
needs to decrease to zero as the rim is approached, we expect that 8h/dR < 0 
for all .R < 0; therefore, the kinetic condition (37) reduces to 


dS ~ ~ 

-^ = Qi{R,t), R<0. 

Furthermore, we have the symmetry conditions 

/9A 

' = 0, Z = 0, R < 0, 


dZ 
90 

^=0, Z = 0, R>0. 
dZ 


(38) 

(39a) 

(39b) 
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By matching to the solutions in region (ii), we obtain the following far-held 
conditions: 


0; = 0s 


6>s(5'(r),0,T). 


(40) 


It is straightforward to see that the bulk equations (35) and the Stefan condition 
(36) are satished by temperatures that are constant in space. Therefore, we have 
that 


0, = 0s= 0 s(^(t),O,t) (41) 

to leading order, which implies the rim moves according to 

f)Q 

— =0s(5(T),r). (42) 

The next-order problem can be used to determine the prohle of the melt near 
the rim; however, this is not required in the subsequent analysis. Finally, we 
note that by matching the melt heights in regions (ii) and (iii), i.e., eh and 
we hnd that 

^(0) ^0, i? ~ Sir). (43) 

We now have all of the ingredients to write down a self-contained problem in 
region (iii). 


3.2.3 A reduced model for 0(1) times 

In region (iii), the melt appears to have zero thickness; it has been collapsed 
onto a circle lying within the Z = 0 plane. The asymptotic matching into the 
inner regions (i) and (ii) provides boundary conditions on this circle. Although 
the melt is effectively treated as having zero thickness, the model still captures 
its evolving shape. 


In region (iii), the temperature field satisfies the equation 


d(h 

dr 


+ 1 , ^ > 0 , 


(44a) 


with 0s = 0 when r = 0. In the far-field, we require that 0s "t as |X| —>■ oo. 
The Z = 0 plane is divided into two regions corresponding to being inside and 
outside of the melt, i? < S(t) and E > S(t), respectively. For points inside of 
the melt, we have Stefan and anisotropic kinetic conditions given by (where we 
drop the (0) subscript on h^^) 


dh 80 s 

dr 



1-k 



Z = 0, R< S{t), 
Z = 0, R< Sir). 


(44b) 

(44c) 
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It should be noted that, because of (44b I, the Stefan condition plays a significant 


role in this regime. This means that the isothermal approximation fails to hold 
and the melt region should no longer be expected to take a kinetic Wulff shape. 
Outside of the melt, we impose a symmetry condition given by 


80 . 


dZ 




= 0 , 


Z = 0, R> S{t). 
The kinetic condition at the rim reads 

Z = 0, R = S{t). 

Finally, it is required that 

HS{t),0,t) = 0 . 


(44d) 


(44e) 


(44f) 


To determine asymptotically consistent initial conditions for the position of the 
rim and the profile of the solid-liquid interface, we examine the early behaviour 


of (441 and match it to the small-time solution given by (24). 


3.2.4 Early behaviour of model for 0(1) times 

The relevant scaling to resolve the early time behaviour and match into the 
small-time regime is given by r = 6. = X = eX, h = eh, and 

S = eS. From the leading-order problem in e, it is straightforward to deduce 
that 9 = T. The leading-order kinetic conditions that hold within the melt and 
at the rim are then given by 


dh 


= r 



8 ^ 

In 


(45a) 

(45b) 


From (24), we see that in the small-time regime, the rim grows like r^/2 -|- 0(e) 
for T ~ 0(1); therefore, we can solve (45b) and by matching we obtain .5(f) = 
f^/2. The solution for S motivates seeking a similarity solution to ( |45a ) of the 
form h = nH{R/n). Using this ansatz in (45a| gives the problem 


HiO-CH'iO = 1 + {H'{0) 


-, 1/2 


(46) 


where C = R/t"^ and H satisfies H{\/2) = 0. The solution is H{C} = ^(1/4 — 
^ 2 ) 1/2 equivalently. 


h = A 



(47) 
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where A = 1 is a constant that can be determined by matching to (24) as 


r ~ 0(1)- From this analysis, we can conclude that the model in (44) should 
have initial conditions for the interface given by 


1 \ 1/2 2 


(48) 


as T ~ 0. For 0 < r <C 1, (48) describes the early growth of the melt in the 


0(1) time regime, which is consistent with the long-term growth in the first 
time regime. 


3.3 Linear Stability for Times of 0(1) 


We now examine the linear stability of the system using the reduced model (44). 


The calculation involves two main steps. First, a base state corresponding to a 
growing axisymmetric melt is computed. Finally, we determine the growth rates 
of small, azimuthally varying perturbations to the base state. Our analysis will 
focus on constructing local solutions valid near, but not too close to, the rim. 

Our calculation of the base state begins by introducing a travelling wave 
coordinate X such that X = R— S{t) and letting Z = Z. We focus on the local 
behaviour of solutions near the rim so that X"^ + Z'^ <C 1. The temperature and 
the melt thickness are written as Og ^ ds{X, Z) and h ~ h{X), where we expect 
from (26) that h{X) ~ hi(—X)^/^ for sufficiently small X. 


Close to the rim, the temperature 9s approximately satisfies Laplace’s equa¬ 
tion: 


d^Os . d‘^9s 


dX^ dZ^ 

The Stefan and kinetic conditions read 
dr dX~ dZ' 


_ - 

dr dX~ " 


= 0 . 



1/2 


Z = Q, X < 0, 
Z = Q, 1<0, 


respectively. The symmetry condition is given by 

^ = 0, Z = 0, X > 0. 

dz 

and the rim evolves according to 

^ = es, z = o, x = o. 

OT 


(49) 

(50) 

(51) 

(52) 

(53) 
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Figure 7: A schematic diagram of the local polar coordinates given hy X = 
u cos (/) and Z = u sin </> that are centred at the rim. It is convenient to write the 
local temperature profile in terms of these coordinates; see text for details. 


Since we have assumed that Og is independent of r, we immediately deduce from 


(53) that the rim moves with a constant velocity, V, given hy V = 0^(0,0). 


An approximate solution for the temperature can be obtained by converting 
to local polar coordinates that are centred at the rim. Thus, we introduce the 
change of variable 

X = ucos(j), Z = usm(j), (54) 

where u is the local radius and (j) is the polar angle measured relative to the 
positive X axis; see Fig. 0(a)- 

An approximate solution for the temperature can be written as 


0s ~ F + cos((j)l2) + yX, 


(55) 


which satishes the symmetry condition ([^ I and where 0i is a constant that can 
be determined from the Stefan condition (50). In particular, by inserting (55) 
in (50) and using the fact that h ^ hi{—X)^>‘^ for A" ~ 0“, we find 




(56) 


so that 01 = Vhi- Using a similar procedure in the kinetic condition (51) shows 


that hi = 1. The parameter 7 is taken to be a free parameter and we will 
investigate the role it plays in controlling the stability of the problem. 

We now investigate the stability of the base state by adding small perturba¬ 
tions of order 5 <C 1 to 0s and S. To simplify matters, we suppose that we are 
looking locally near {X, Y, Z) = (Ur + X, U, Z), where X"^ + Y'^ + Z'^ and 
can consider the rim as a straight line on these scales. Taking the rim to be flat 
is reasonable when the perturbation wavenumber in the azimuthal direction is 
large. Note that X = Y = Z = 0 corresponds to a point on the base-state rim 
and, thus, we have effectively attached a Cartesian coordinate system to this 
point. We write the local temperature and the position of the rim as 


0s ^ 0s(X, Z) + (50s(l, Z) 
S'~ Ur + 


(57a) 

(57b) 
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where k and m denote the wavenumber and growth rate of the perturbations, 
respectively, and 9s is given by (551. The perturbation to the temperature 
satishes the equation 


d^Qs 

dX^ 


d'^Qs 

dZ^ 


- k^Qs = 0 , 


together with 


^ = 0, Z = 0, X >0. 
dz 


(58) 


(59) 


The solution can be found using the local polar coordinates in (541 and written 
as 


0s = 0iu ^/^e cos(((>/2), 


(60) 


where 0i is a constant that is to be determined. An equation governing the 


perturbation to the rim position can be derived from inserting (571 into the 
kinetic condition dS/dr = 9 s(S{t),t), expanding about d <C 1, and taking the 
0{6) part: 


We note that 


f)f) 

mS=^S + &s, Z = 0,X = 0. 
dX 


oX ^ 


(61) 


(62) 


as X ~ 0 and Z = 0, both of which become singular as X —)■ 0. In order for the 


kinetic condition (61) to remain well defined, we need 0i = —{V/2)S^ which 
yields 

mS = jS, (63) 

i.e., the perturbation growth rate m is exactly equal to the parameter 7 in the 


base-state temperature profile (55). This linear analysis thus indicates instabil¬ 


ity if there is a background temperature gradient in the direction of propagation 
of the rim, 7 > 0, but stability for a negative gradient, 7 < 0. Note that in the 
case of instability, the growth rate of the perturbations is independent of the 
wave number, in contrast to unstable Hele-Shaw or Stefan problems without 
surface tension/energy, where growth rate increases with wave number and can 
be arbitrarily high. Note that similar stability results for another free boundary 
problem were obtained in Howison et al. [Si- 

Given the absence of exact and of approximate long-time solutions about 
which to perturb, it is not immediately apparent what values 7 might take in 
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practice. Intuitively we might expect 7 to be positive, since melting at the 
interface has the effect of locally reducing temperature, at least for relatively 
low times r. The simulations by Hennessy [15] support this claim, although 
they do not consider heat transfer in the axial direction. If 7 is positive, we 
then expect a mild instability whose form will also be influenced by any further 
anisotropy, for instance, the usual six-fold one in the {X, Y) plane. 


3.4 Other Anisotropy Functions 

We now briefly outline the results that are obtained for the anisotropy functions 
Q and Full details about the solutions in the early-time regime and the 
solution of Charpit’s equations are given in Appendix [^ 

3.4.1 Dynamics with /{ip) = e + sin^ ip 

In the early-time regime given by r = the solid-liquid interface can be 

written parametrically as 

X = [Q;-|-s(I-|-e-|- sin^ </?)] cos tp, ^ = [a -I- s(e — cos^ (/?)] sin tp, (64) 

where (p G [0, 27r) and s = r^/2. Interface profiles at various times are shown in 
Fig. I (a). The interface remains smooth until s = r^/2 = a/(I — e), at which 
point a corner develops at the rim due to intersecting characteristics. The early 
growth of the rim for s/a e scales like S' ~ r; however, the longer-term growth 
of the rim for s/a ^ e is reduced by the corner and we find that S ~ 

The thickness of the melt in the axial direction grows like er^ for all times. 


For larger times, the separation of length scales in the radial and axial di¬ 
rections can, in principle, be exploited and the model can be reduced using a 
similar analysis to that in Sec. |3.2| However, the current model is expected to 
require additional mechanisms such as surface energy to act to regularise the 
corner. Therefore, we do not proceed with the model reduction in this case. 
Nevertheless, we note that because of the slower radial growth for early times, 
in getting to terms to balance in a model equivalent to (44), larger time and 
temperature scalings are needed: r = and 0 = . 


3.4.2 Dynamics with /('(/') = e/(l + e ~ sin^ ip) 

In this case, the position of the interface in the early-time regime, t = 0{a^^^), 
is given by 


A = 


a -I- 


es(3sin^ (fi + e) 
(sin^ p + 


cos tp. 


Z = 


es(3cos^ p — 1 — e) 


(sin^ p ■ 


sin(/j, (65) 


where, again, s = r^/2. Fig. [^(b) shows the corresponding interface profiles at 
various times. Here, the corner appears in the very early stages of melt growth, 
in particular, when s = r^/2 = eo/(2 — e). The growth of the rim scales like 
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Figure 8: Early-time evolutions of an initially spherical melt for anisotropy 
functions /(V') = e + sin^ tp (panel a) and /{tp) = e/(l + e — sin^ ip) (panel 
b) when e = 0.1. The curves in panels (a) and (b) are obtained from the 
solutions ( |64| and (65), respectively. The positions of the interface are shown 
at equally spaced values of s given by s = 0, 2, 4, 6, 8, and 10, corresponding 
to rescaled dimensionless times t = 0, 2, 2.83, 3.46, 4.0 and 4.47, respectively. 
Both anisotropy functions lead to the formation of a corner, and this happens 
when f = 1.49 in panel (a) and f = 0.32 in panel (b). As f —)■ oo, the interface 
profiles approach the kinetic Wulff shapes shown in Fig. 
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S' ^ a + (ae)^/^T for s a/e and like S ^ er^ for s a/e. The axial growth 
scales like er^ for all time. 

For this particular anisotropy, the eventual growth of the melt both parallel 
and normal to the c axis is the same order of magnitude. The aspect ratio of 
the melt roughly approaches 5:2 and, therefore, it is not possible to simplify the 
model for 0(1) times. 


3.4.3 Commonalities of the early-time growth 

The three anisotropy functions that we consider produce interface prohles with 
common growth features in the early-time regime. For instance, all three cases 
lead to melts that evolve into their kinetic Wulff shapes given by ([^. In fact, 
an analysis for arbitrary anisotropy functions in Appendix shows this will 
always be the case. Furthermore, the growth of the melt in the axial direction, 
i.e. along the c axis, is always found to be quadratic with time. As shown in 
Appendix]^ if nucleation occurs much later than when the system is irradiated, 
then the axial growth becomes linear for all anisotropy functions. 

4 Discussion and Conclusion 

In this paper, we have formulated and analysed a mathematical model describ¬ 
ing the anisotropic growth of a Tyndall figure into a crystal of superheated ice. 
Both the solid and liquid phases are assumed to be volumetrically heated by 
the absorption of incoming radiation, which drives the melting process. The 
anisotropic growth of the Tyndall figure is a result of the molecularly smoothly 
basal planes of the ice crystal melting at a much slower rate in comparison to 
molecularly rough prism planes. This phenomenon is modelled using a kinetic 
coefficient that depends on the orientation of the solid-liquid interface. The 
relationship between the kinetic coefficient and the crystal orientation is quan¬ 
tified through an anisotropy function. Our analysis indicates that there are two 
key time regimes for the melt evolution. The first of these describes the rapid 
initial growth of the Tyndall figure into its kinetic Wulff shape due to volumet¬ 
ric heating. The second time regime describes the slower, diffusion-dominated 
growth. 

The problem in the first time regime amounted to solving an anisotropic 
Eikonal equation. Remarkably, it was possible to obtain an analytical solu¬ 
tion to this equation for an arbitrary anisotropy functions. Using this solution, 
we examined the interface profiles and kinetic Wulff shapes that are obtained 
for three different anisotropy functions. These anisotropy functions led to a 
rich variety of melt shapes including long rectangles with rounded ends, oblate 
spheroids, as well as thick and thin lenses. Qualitatively, we found that the 
smoothest melts and the smallest aspect ratios occur when the anisotropy func- 
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tion has broad maxima; anisotropy functions with narrow maxima gave rise to 
corners and lens-shaped melts that can have order-one aspect ratios. Regard¬ 
less of the anisotropy function, the thickness of the melt in the direction of the 
c axis was found to grow quadratically with time. This is in contrast to the 
radial growth parallel to the basal planes, which was highly dependent on the 
anisotropy function. These findings have important practical implications, as 
they suggest that experimental data for the radial growth of the melt can aid 
in the determination of likely candidates for the anisotropy function. This is 
not the case for axial growth, which is predicted to be roughly the same for all 
anisotropy functions. 


By exploiting the thin aspect ratio of the melt figure, we showed that a 
simplified model for the evolution in the second time regime can be derived by 
systematically collapsing the three-dimensional melt figure to a two-dimensional 
surface with zero thickness along the axial direction. This model was then used 
to carry out a linear stability analysis, the results of which suggest that an 
instability will occur if the temperature field locally increases in the direction of 
radial growth. Such an instability would likely lead to fingers and could drive 
the formation of a Tyndall star similar to that shown in Fig. [l] 


The results from our analysis, in combination with the experimental ob¬ 
servations by Takeya [29j . may give some insight into appropriate anisotropy 
functions for the melting of ice crystals. In particular, the melts documented 
by Takeya have a remarkably constant thickness in direction of the c axis which 
diminishes relatively rapidly near the rim. In addition, the aspect ratio of the 
melt is small and on the order of 1:10. These observations suggests that an ap¬ 
propriate anisotropy function for modelling the growth of Tyndall figures would 
be similar to that in ([^ but with much broader maxima at '0 = ±7r/2. 


Further predictions about Tyndall star evolution can be accessed through 
numerical simulations of our model. From a computational perspective, simpli¬ 
fied models such as (44) are advantageous due to the reduced dimensionality of 
the free boundary and are relatively straightforward to implement. Numerical 
simulations of such a model can provide insights into when the condition for in¬ 
stability is satisfied and offer a means of probing nonlinear melt morphologies. 
Furthermore, such simulations could explore whether the onset of instability is 
linked to growth along the c axis, which has been suggested by experimental 
studies [IHlIini- Thus, there is a wide range of exciting and unanswered prob¬ 
lems relating to the formation and evolution of Tyndall stars, and we hope this 
work not only provides some of the foundations that can aid in tackling these, 
but also motivation for doing so. 
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A Solution of Charpit’s Equations for the Anisotropic 
Eikonal Equations 

An asymptotic analysis of the model revealed that the early-time interface pro¬ 
files can be obtained by solving an anisotropic Eikonal equation of the form 


(s|+ (sx(sx+ 4) (66a) 


where sx = ds/dX and sz 
condition 


dsjdZ. Equation (661 is supplemented with the 


So — s(Ao, Zq) = 0, Xq + Zq — c?. 


(66b) 


The solution to this problem can be obtained using Charpit’s equations, which 
generalise the well-known method of characteristics to fully nonlinear first-order 
hyperbolic partial differential equations (PDEs) [57]. We recall that when ap¬ 
plying the method of characteristics, one must simultaneously solve for the 
characteristic directions and the solution to the PDE on these characteristics. 
The idea behind Charpit’s method is to treat the first derivatives of the solution 
to the PDE as additional unknowns that must be found along the characteristic 
directions. Thus, when applying Charpit’s method to this problem, we must 
simultaneously solve for the characteristic directions, X and Z, as well as the 
solution a and its derivatives sx and az along the characteristics. Although 
these five unknowns are effectively treated as independent variables, Charpit’s 
equations ensure that they always vary in a consistent manner. 


the PDE in (66a) as 


To apply Charpit’s method to (66), we first let p = ax, q = sz, and we write 


G{X, Z, a,p, q) = (p2 + g2)l/2; ^ q 2 yl/ 2 ^ _ ^ ^ 


(67) 


The condition in (66b I can be treated as initial data and parametrised according 
to 


■so(<7’) = s{Xoi<p),Zo{(p)) 
Xo{ip) = a cos ip, 

Zo{p) = asintp. 


0, c = o, 

(68a) 

C = o, 

(68b) 

C = o, 

(68c) 
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where ip € [0, 2tt) and C is an arbitrary parameter that measures distance along 
each characteristic direction. Initial conditions for p and q, given by po and go, 
can be obtained by (i) differentiating the condition so{ip) = s{Xq{(p), Zo(ip)) = 0 
with respect to ip and (ii) requiring the PDE (67) to hold on the initial curve, 
G{Xq, Zq, sqjPq, go) = 0. By simultaneously solving two these equations, we 
obtain 


Poi^p) 


COS ip 

/(cos ip) ’ 


qo{p) 


sin (f 
/(cos ip) ’ 


C = 0. 


(68d) 


Charpit’s equations for this problem can be written as 


. dG 
dp' 

(69a) 

dG 



(69b) 

s = 1, 

(69c) 

0 = 0, 

(69d) 

9 = 0, 

(69e) 


where the dot denotes differentiation with respect to Upon solving these 
equations with the initial conditions in (68), we find that s = ^, so that ^ can 
be replaced by s. In addition, we have p = pq, q = q^, and 


X = [a-\- s/(cos ip)] cos ip + s/'(cos p) sin^ p, (70a) 

Z = [a + s(/(cos p) — /'(cos (/?))] sin p, (fOb) 


with the prime denoting derivative with respect to argument. 


For the anisotropy function (a) fipp) = (e^ + sin^ we have that f{w) = 

(e^ + After inserting this expression into ( [70| and some algebra, the 

solution can be written as 


a: = 



^(1 + e ^) \ 

(e^ + cos^ p)^i'^ ) 


cos p, 


Z = [ a + 


cos- 




sm p. 


(71) 


The properties of this solution are described in Sec. |3.1| For the anisotropy 
functions (b) f{ip) = e + sin^ 0 and (c) /(0) = e/(l + e — sin^ 0), we find that 


X = [a + s(l + sin^ p + e)] cos p, Z = [a — s(cos^ p — e)] sin p, (72) 


and 

a: = 


es(3sin^ p + t) 


cosp^ Z = 


es(3cos^ p — 1 — e) 


(sin^ p + e)"^ j ’ [ (sin^ p + e)^ 

respectively. These solutions with e = 0.1 are shown in Figs. and 


sin(/j, (73) 
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For case (b), focusing on that part of the free boundary lying in the first 
quadrant, 0 < (/? < 7r/2, we see that some of the characteristics are directed 
down, towards Z = 0, and intersection of characteristics starts, on the X axis, 
when s = a/(l — e)^a at X = a + a(l + e)/(l — e) = 2 q;/( 1 — e) ~ 2a. For 
later times this method of characteristics indicates multiple-valued solutions. 
To avoid this, the convex part of the curve is taken, giving corners on Z = 0 for 
s > a/(l — e). These would be expected to be rounded off by any sort of surface- 
tension or surface-energy effects so that a Gibbs-Thomson term is introduced 
into the free-boundary conditions. A mathematically simpler way of regularising 
the problem would be to replace the anisotropic Eikonal equation, which is a 
first-order flow, by a mean-curvature flow. Results of Barles & Souganidis [3] 
could be applied to give continuous dependence of solutions on the coefficient of 
any curvature term included in (14). This would again indicate that we should 
get the interface by taking the convex part of the curve. 


The same corner formation is seen for the anisotropy function (c). In this 
case, the corner forms very quickly, when s = ea/(2 — e), and close to the initial 
free boundary, at X = 20/(2 — e). 


The range of possible short-time interface behaviour is large because the 
growth in the X direction can have quite different qualitative behaviour. For 
case (b), with s <C a/e, X = 2s(e + ~ 2(as)^/^ = (2a)^/^T so that the 

growth is only linear in time. The final case (c), has cos^ v? ^ 1 — (2es/a)^/^ 
and X ^ q;(1 -|- (2es/a)^/^) = a(l -I- (e/af-l'^T). 

For large times, in the sense of s a/e, the behaviour of Z is the same for 
all three anisotropy functions: Z es ^ eT^ 12. However, the long-time growth 
in the X direction is reduced, thanks to the appearance of the corner. For (b), 
the corner’s position is, in general, given by X = [a-|-s(l-|-sin^ tp + e)] cosp with 
Z = [a — s(cos^ (f — e)] sin(/j = 0. Since 0 < (p < 7r/2, this gives cos^ p = e + a/s 
and X = 2s(e-|-a/s)^/^ ^ 2e^/^s = e^l'^T^ for s 3> a/e. Very similar calculations 
for (c) show that the corner location can be obtained implicitly from 

X 2cos(/ 3 es (1 — cos^ 
a 3 cos^ Lp a 3 cos^ ip — 1 

for 0 < (^ < cos“^(l/-\/3). For s S> a/e, this gives cos^ ~ (1/3) and we get 
X ~ 33/V2es = (33/74)eT2. 


By taking the modified time variable s sufficiently large in comparison to a 
in (70), the longer-term interface profile for an arbitrary anisotropy function is 
given by 


X/s ~ /(cos p) cos p + //cos p) sin^ p, (75a) 

Z/s ~ /(cos p) sin(^ — /'(cos p) cosy>siny), (75b) 
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independent of the initial shape. Equation (75) is, in fact, equivalent to ([^ and 
therefore, the interface profiles approach the kinetic Wulff shapes. The direction 
of the characteristics, ZjX^ can be differentiated with respect to Lp to check if 
this ever decreases, leading to corner formation from an initially convex shape. 
After some manipulation, the derivative turns out to be 




The criterion for a continued smooth interface is then / + d^//d'0^ > 0. Cahoon 
et al. [ 6 ] and Wettlaufer et al. [34j find the same basic law for interface motion 
gives the rate of change dre/ds = (/+ d^//di/'^)K^ for the interface curvature k. 
The same key combination appears in curvature-flow models for phase change 
with significant Gibbs-Thomson effect mna. In these works, the {f+d'^f/d-tp'^) 
term multiplies curvature in the velocity law and, to avoid negative diffusion, 
all angles making (/ -|- d'^f/dtp'^) positive are prohibited, leading to corners in 
the interface for all positive time. 


B The Role of Nucleation 


We now give a brief discussion of the effect of surface energy in the nucleation 
process, while still neglecting the air bubble that appears in the melt. We con¬ 
centrate on the implications of the balance between the superheat temperature 
and the local equilibrium temperature for a spherical liquid body of a given size; 
Chadham et al. [7] discusses related effects in the growth of crystals when the 
Gibbs-Thomson effect is the only stabilising action. 

We suppose that nucleation occurs when the temperature in the solid exceeds 
a nucleation temperature T„ given by the Gibbs-Thomson relation 


— Tn 



(76) 


where 7 is the interface energy, a„ is the nucleation radius. The time at which 
nucleation occurs, measured relative to the moment the system is irradiated, is 
denoted by t„. Before nucleation occurs, the temperature in the solid increases 
like Ts = Tq + qgt/{pcps)] therefore, the nucleation time and temperatures can 
be related via = pCps{Tn — TQ)/qs. 

So far we have been assuming that the nucleation temperature is close to the 
bulk melting temperature, Tn — Tq, so that nucleation immediately occurs upon 
irradiation, resulting in an initial liquid-solid interface that is approximately a 
sphere of radius a, which is small compared to e£. The condition a <C allows 
the melt to become a developed spheroid when the dimensionless time r is 0 ( 1 ) 
in size. Note that if e <C a ^ 1 there is a significant change to Sec. |3.2[ with 
the Tyndall figure no longer being of thickness order e. 

We now consider the opposite case whereby T„ Tg so that nucleation 
occurs much later than when the system is irradiated. The bulk temperatures in 
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this case will be large during the early evolution of the melt and will influence its 
growth kinetics. To study the behaviour in this late-nucleation regime, we non- 
dimensionalise by writing x = IX, t = t„ + (P/Ks)f, and T = Tq + XTO, 

where AT = T„ — Tq = 2"f/{pLan) and I = ks/{pCpsATK). The dimensionless 
volumetric heat sources given by qii'^ / {pCpsAT) characterise the temperature 
rises that occur on the diffusive time scale due to absorption relative to the 
nucleation temperature AT = T„ — Tq. These relative temperature rises are 
expected to be small so the volumetric source terms are neglected from the 
model, i.e., we take qiP/{pCpgAT) ~ 0. The dimensionless bulk equations for 
the temperatures can be written as 

- - - 

^ X e (77a) 

BO 

Cp-^ = kV%, (77b) 

which have initial conditions 9s = 0i = 1 when r = 0 and far-held conditions 
ds ~ 1 for |X| —)■ c». At the free boundary, the Stefan condition reads 

where the Stefan number is now given by /3 = L/[cpgAT). The anisotropic 
kinetic condition is 


v = 9jf(^), Xef(f), (77d) 

The initial interface r(0) is assumed to be a circle of dimensionless radius a = 
anil. 

In order to obtain the same asymptotic regimes as in the early-nucleation 
case considered in the main text, we let j5~^ = be and require the dimensionless 
initial melt radius to satisfy d e. The condition Tn S> Tg imposes an addi¬ 
tional restriction on the dimensionless nucleation radius given by d <C 2'^/ [pLtj. 
Thus, in dimensional terms, we require 


a„ <C min 




(78) 


We now summarise the early-time, r <C d, and order-one time, r = 0(1), 
problems in the late-nucleating regime. 

The early-time problem valid for f <C d can be deduced by repeating the 
analysis of Sec. |3.1[ The lack of a volumetric heat source means that the 
leading-order temperatures (in d) are constant in time, 9i = Og = \. Thus, 
the anisotropic kinetic condition becomes v = fipp), which is now autonomous 
in the time variable r. As a consequence, the growth kinetics of the melt are 
modified. For the anisotropy function given by /(^) = (e^ -I- sin^ '0)^^^) we find 
that the interface can be written parametrically as 


A = 


f(l + a 


(e^ -I- cos^ (/5 )i/2 


cosp, Z = 


(e^ -I- cos^ "PpA 


sin(^, (79) 
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where (p G [0, 27r). Thus, the thickness and rim of the melt now grow linearly 
with time rather than quadratically. However, the morphological characteristics 
of the interface remain the same as in the early-nucleation regime and, in par¬ 
ticular, the kinetic Wulff shapes are still approached in the longer term. Similar 
changes are seen for other anisotropy functions as well; that is, the powers of r 
in the growth laws are reduced by a factor of two. 

For 0(1) times and the anisotropy function / = (e^ -|-sin^ tti® analysis 

in Sec |3.2| can also be repeated in order to derive a simplified model that collapses 
the melt region onto the Z = 0 axis. In particular, the temperature in the solid 
satishes the equation 


Z>0, (80a) 

with dg = 0 when f = 0 and 0s ~ 1 as |X| —>■ oo. The boundary conditions on 
Z = 0 are given by 


dh - 89s 



1-h 


55 

df 


= 9 . 


dZ 


= 0 , 



Z = 0, R< 5(f), 

Z = 0, R< 5(f), 

z = o, 5 = 5(f), 
Z = 0, R> 5(f). 


Finally, we require that 


h{S{T),0,T) = 0. 


(80b) 

(80c) 

(80d) 

(80e) 

(80f) 
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